home *** CD-ROM | disk | FTP | other *** search
/ Compendium Deluxe 2 / LSD and 17bit Compendium Deluxe - Volume II.iso / a / prog / asmsrc / thesource-7.lha / Source / NumericalMethods.lha / NumericalMethods / sqrt.c < prev    next >
C/C++ Source or Header  |  1994-05-27  |  9KB  |  382 lines
  1. /* --------------------------------- sqrt.c --------------------------------- */
  2.  
  3. /* This is a collection of sqrt(ulong) programs off the net (see later for the
  4.  * original notes). I added my own version, eyal(), which uses the standard
  5.  * iterative approach. eyal0() is the simple one while eyal() uses a table
  6.  * lookup for the initial guess.
  7.  *
  8.  * Some functions were adjusted so that all truncate in the same way.
  9.  *
  10.  * Eyal Lebedinsky    (eyal@ise.canberra.edu.au)
  11. */
  12.  
  13. #if 0
  14.  
  15.  These are the results for running the program on a 486DX50 on
  16. msdos(quickC), msdos(C7.00) and Linux(gcc2.2.2). Note that the ftime on
  17. Linux gives 10ms resolution. The table shows time in 1ms units. All
  18. compiles were done with maximum optimization (best speed options).
  19.  
  20.  The summary line shows what the compiler can contribute (compare qc and c7)
  21. and what the machine capability can do (gcc test done in 32-bit mode).
  22.  
  23.  
  24. calling sqrt() 100000 times:
  25.  
  26. function    gcc     c7      qc     check
  27.                                              
  28. null         60     84     120         0 
  29. null         60     84     120         0 (acid)
  30. rodent      350  2,818   3,165  21032170 
  31. rodent      350  2,921   3,275  74047927 (acid)
  32. grupe       690  3,364   3,681  21032170 
  33. grupe       700  3,516   3,826  74047927 (acid)
  34. dj          980  4,389  19,681  21032170 
  35. dj        1,050  4,582  20,239  74047927 (acid)
  36. thyssen   1,100  4,462  20,282  21032170 
  37. thyssen   1,160  4,643  20,671  74047927 (acid)
  38. kskelm    3,030 14,878  14,625  21032170 
  39. eyal0       700  2,099   2,464  21032170 
  40. eyal0       830  2,696   3,134  74047927 (acid)
  41. eyal        320    555     802  21032170 
  42. eyal        350    625   1,405  74047927 (acid)
  43.          ====== ====== =======
  44.          11,730 51,716 117,490
  45.  
  46.  
  47. Original posting:
  48. =================
  49.          
  50. $From comp.sys.ibm.pc.programmer Tue Oct  8 09:16:35 1991
  51. $Newsgroup: comp.sys.ibm.pc.programmer/3360
  52. $Subject: Summary: SQRT(int) algorithm (with profiling)
  53. $From: warwick@cs.uq.oz.au (Warwick Allison)
  54. $Sender: news@cs.uq.oz.au
  55. $Date: 7 Oct 91 01:20:38 GMT
  56. $Message-Id: <4193@uqcspe.cs.uq.oz.au>
  57. $Path: csc.canberra.edu.au!manuel!munnari.oz.au!bunyip.cc.uq.oz.au!uqcspe!cs.uq.oz.au!warwick
  58. $Reply-To: warwick@cs.uq.oz.au
  59. $Followup-To: comp.sys.amiga.programmer
  60. $Lines: 177
  61.  
  62.  
  63. Thanks to all those who responded.
  64.  
  65. Profiles:
  66.  
  67.  %time  cumsecs  #call  ms/call  name
  68.   52.4     3.16  99999     0.03  _kskelm
  69.   12.1     3.89  99999     0.01  _thyssen
  70.   10.8     4.54  99999     0.01  _grupe
  71.   10.6     5.18  99999     0.01  _dj
  72.    5.3     5.98  99999     0.00  _rodent
  73.  
  74. All gave the same accuracy except for grupe, which also rounded rather
  75. than truncating.
  76.  
  77. (I added this to rodent() at slight cost)
  78.  
  79. Thanks again all,
  80. Warwick
  81. --
  82.   _-_|\       warwick@cs.uq.oz.au
  83.  /     *  <-- Computer Science Department,
  84.  \_.-._/      University of Queensland,
  85.       v       Brisbane, AUSTRALIA.
  86.  
  87. ---------------- code below -----------------
  88. #endif
  89.  
  90. #include <stdio.h>
  91. #include <stdlib.h>
  92.  
  93. #if 1
  94. /* This is for msdos
  95. */
  96. #include "tick.h"
  97. #define GET_TIME    (GetLowTickCount()/10L)
  98. #endif
  99.  
  100. #if 0
  101. /* This is for unix
  102. */
  103. #include <time.h>
  104. #include <sys/timeb.h>
  105.  
  106. #define GET_TIME    timer_milli()
  107.  
  108. static unsigned long
  109. timer_milli (void)
  110. {
  111.     struct timeb    tm;
  112.  
  113.     ftime (&tm);
  114.     return (tm.time*1000L + tm.millitm);
  115. }
  116. #endif
  117.  
  118. /*
  119.     This integer sqrt function is about *15* times faster than
  120.     the standard double-precision math function under lattice.
  121.     It was the result of a thread on comp.graphics a while ago.
  122.     I've tested the precise accuracy of it's results up to 600000.
  123. */
  124.  
  125. static unsigned long rodent(v)
  126. unsigned long v;
  127. {
  128.     register long t = 1L<<30, r = 0, s;
  129.  
  130. #define STEP(k) s = t + r; r >>= 1; if (s <= v) { v -= s; r |= t;}
  131.  
  132.     STEP(15);   t >>= 2;
  133.     STEP(14);   t >>= 2;
  134.     STEP(13);   t >>= 2;
  135.     STEP(12);   t >>= 2;
  136.     STEP(11);   t >>= 2;
  137.     STEP(10);   t >>= 2;
  138.     STEP(9);    t >>= 2;
  139.     STEP(8);    t >>= 2;
  140.     STEP(7);    t >>= 2;
  141.     STEP(6);    t >>= 2;
  142.     STEP(5);    t >>= 2;
  143.     STEP(4);    t >>= 2;
  144.     STEP(3);    t >>= 2;
  145.     STEP(2);    t >>= 2;
  146.     STEP(1);    t >>= 2;
  147.     STEP(0);
  148.  
  149. /*  if (r<v) return r+1;    add for rounding */
  150.  
  151.     return r;
  152. }
  153.  
  154.  
  155. static unsigned long grupe(x)
  156. unsigned long x;
  157. {
  158.         register unsigned long xr;      /** result register **/
  159.         register unsigned long q2;      /** scan-bit register **/
  160.         register int f;                 /** flag (one bit) **/
  161.  
  162.         xr = 0;                         /** clear result **/
  163.         q2 = 0x40000000L;               /** higest possible result bit **/
  164.         do
  165.         {
  166.           if((xr + q2) <= x)
  167.           {
  168.             x -= xr + q2;
  169.             f = 1;                      /** set flag **/
  170.           }
  171.           else f = 0;                   /** clear flag **/
  172.           xr >>= 1;
  173.           if(f) xr += q2;               /** test flag **/
  174.         } while(q2 >>= 2);              /** shift twice **/
  175. /*      if(xr < x) return xr +1;         add for rounding */
  176.         return xr;
  177. }
  178.  
  179.  
  180. static unsigned long kskelm(val)
  181. unsigned long val;
  182. {
  183.         register unsigned long rt = 0;
  184.         register unsigned long odd = 1;
  185.  
  186.         while (val >= odd) {
  187.          val -= odd;
  188.          odd += 2;
  189.          rt += 1;
  190.         }
  191.         /* sqrt now contains the square root */
  192.         /* val now contains the remainder */
  193.     return rt;
  194. }
  195.  
  196. static unsigned long dj(val)
  197. unsigned long val;
  198. {
  199.   unsigned long result = 0;
  200.   unsigned long side = 0;
  201.   unsigned long left = 0;
  202.   int digit = 0;
  203.   int i;
  204.   for (i=0; i<16; i++)
  205.   {
  206.     left = (left << 2) + (val >> 30);
  207.     val <<= 2;
  208.     if (left >= (side<<1) + 1)
  209.     {
  210.       left -= (side<<1)+1;
  211.       side = (side+1)<<1;
  212.       result <<= 1;
  213.       result |= 1;
  214.     }
  215.     else
  216.     {
  217.       side += side;
  218.       result <<= 1;
  219.     }
  220.   }
  221.   return result;
  222. }
  223.  
  224. static unsigned long thyssen(val)
  225. unsigned long val;
  226. {
  227.   register unsigned long result = 0;
  228.   register unsigned long side = 0;
  229.   register unsigned long left = 0;
  230.   int i;
  231.  
  232.   for (i=0; i<sizeof(unsigned long)*4; i++) /* once for every other bit */
  233.   {
  234.     left = (left << 2) + (val >> 30);
  235.     val <<= 2;
  236.     if (left >= side*2 + 1)
  237.     {
  238.       left -= side*2+1;
  239.       side = (side+1)*2;
  240.       result <<= 1;
  241.       result |= 1;
  242.     }
  243.     else
  244.     {
  245.       side *= 2;
  246.       result <<= 1;
  247.     }
  248.   }
  249.   return result;
  250. }
  251.  
  252. static unsigned long eyal0(x)        /* used for initialization only */
  253. unsigned long x;
  254. {
  255.     register unsigned long    r, t;
  256.     register long        e;
  257.  
  258.     if (x & 0xffff0000L)
  259.         r = 662 + x / 17916;
  260.     else if (x & 0x0000ff00L)
  261.         r = 3 + x / 70;
  262.     else
  263.         r = 2 + x / 11;
  264.  
  265.     do {
  266.         t = x / r;
  267.         e = (long)(r - t) / 2;
  268.         r = (r + t) >> 1;
  269.     } while (e);
  270.     return (r);
  271. }
  272.  
  273. static unsigned int    sqrtab[256+1] = {0};
  274.  
  275. static void
  276. set_eyal ()
  277. {
  278.     int    i;
  279.  
  280.     for (i = 0; i < 256; ++i)
  281.         sqrtab[i] = eyal0 (i*256L*256L*256L);
  282.     sqrtab[256] = 0xffffU;
  283. }
  284.  
  285. static unsigned long eyal(x)
  286. unsigned long x;
  287. {
  288.     register unsigned int    r, t;
  289.     register int        e;
  290.  
  291. /* Select the intial guess. Ensure that it is ABOVE the qsrt so that the
  292.  * long/short divide will fit into a short (or else we get a divide
  293.  * overflow).
  294. */
  295.     if (x >= 0x00010000UL)
  296.         if (x >= 0x01000000UL)
  297.             if (x >= 0xfffe0001UL)
  298.                 return (0xffffU);
  299.             else
  300.                 r = sqrtab[(x>>24)+1];
  301.         else
  302.             r = sqrtab[(x>>16)+1] >> 4;
  303.     else if ((unsigned int)x >= 0x0100U)
  304.         r = sqrtab[((unsigned int)x>>8)+1] >> 8;
  305.     else
  306.         return (sqrtab[x] >> 12);
  307.  
  308. /* Iterate until error is zero.
  309.  * It was measured that on randomly large numbers (acid test) we go round
  310.  * the loop on average 2.2 times.
  311. */
  312.     do {
  313.         t = (unsigned int)(x / r);
  314.         e = (int)(r - t) >> 1;
  315.         r -= e;
  316.     } while (e);
  317.  
  318.     return (t);
  319. }
  320.  
  321. static unsigned long null(x)
  322. unsigned long x;
  323. {
  324.     return (x);
  325. }
  326.  
  327. static void
  328. test (func, name, n, acid)
  329. unsigned long (*func) ();
  330. char        *name;
  331. unsigned long    n;
  332. int        acid;
  333. {
  334.     unsigned long    i, x, xx, t, a;
  335.  
  336.     a = 0L;
  337.     x = acid ? 0xffffffffUL/n : 1L;
  338.     t = GET_TIME;
  339.     for (xx = x, i = 0L; i < n; i++, xx += x)
  340.         a += (*func) (xx);
  341.     t = GET_TIME - t;
  342.     printf ("%-10s %7ld %10ld %s\n", name, t, a, acid ? "(acid)" : "");
  343.     fflush (stdout);
  344. }
  345.  
  346. #define TEST(f,n,acid) \
  347.     test (f, n, loops, 0);    \
  348.     if (acid) test (f, n, loops, 1)
  349.  
  350. int
  351. main (argc, argv)
  352. int    argc;
  353. char    *argv[];
  354. {
  355.     long    i, a, t, loops;
  356.  
  357.     if (argc) {
  358.         loops = atol (argv[1]);
  359.         if (loops <= 0) {
  360.             printf ("bad count\n");
  361.             exit (1);
  362.         }
  363.     } else
  364.         loops = 10000L;
  365.  
  366.     printf ("calling sqrt() %ld times:\n", loops);
  367.  
  368.     printf ("\n%-10s %7s %10s\n\n", "function",  "time", "check");
  369.  
  370.     TEST(null,"null",1);
  371.     TEST(rodent,"rodent",1);
  372.     TEST(grupe,"grupe",1);
  373.     TEST(dj,"dj",1);
  374.     TEST(thyssen,"thyssen",1);
  375.     TEST(kskelm,"kskelm",0);    /* too slow for acid test */
  376.     TEST(eyal0, "eyal0",1);
  377.     set_eyal ();
  378.     TEST(eyal,"eyal",1);
  379.  
  380.     exit (0);
  381. }
  382.